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Abstract 

The objective of this paper is to describe a new method for identification of a 
continuous-time multi- input and multi-output bilinear system. The approach is to make 
judicious use of the linear- model properties of the bilinear system when subjected to a 
constant input. Two steps are required in the identification process. The first step is to 
use a set of pulse responses resulting from a constant input of one sample period to 
identify the state matrix, the output matrix, and the direct transmission matrix. The 
second step is to use another set of pulse responses with the same constant input over 
multiple sample periods to identify the input matrix and the coefficient matrices 
associated with the coupling terms between the state and the inputs. Numerical examples 
are given to illustrate the concept and the computational algorithm for the identification 
method. 


Introduction 

System identification is a methodology used to characterize a dynamical or other 
engineering system with measurements of the input-output signals. Mathematicians and 
engineers have developed a number of approaches to address the identification problem. 
The identification of a linear time- invariant system is relatively veil understood and 
theoretically well developed [1,2]. This is not true for the identification of a nonlinear 
system, although some progress has been made in the identification of nonlinear systems 
over the past few decades [3-19], 

There is a class of nonlinear systems called bilinear systems whose dynamics are 
jointly linear in the state and the force variables. It is a simple nonlinear extension of a 


linear system. The concept of bilinear systems was introduced in the 1960’s (see the 
surveys of Refs. [5] and [6]). References [7] and [8] provide a survey of bilinear-related 
system- theory methods and their contributions to problems such as stabilization, 
controllability, and observability. Bilinear systems have been studied extensively and 
applied successfully to several problems [15]. Recently, research activities in 
identification of bilinear systems have been focused on the so-called “discrete-time” 
model identification [1 9]. The discrete-time model is an approximation obtained by 
linearizing the continuous one with a method such as the finite difference. In contrast, we 
focus on the identification of a continuous-time bilinear system without any 
approximation. 

A new method is introduced in this paper for identification of a continuous-time 
multi- input and multi-output bilinear system. When the input of a bilinear system is a 
constant, the bilinear system becomes a linear system. This special characteristic is the 
basis for the identification method. Two steps are required for the identification process. 
The first step begins with generating a set of pulse responses with a constant input 
applied one at a time over one sample period. The pulse responses are then used to form 
a Hankel matrix consisting of system Markov parameters to identify the state matrix, the 
output matrix, and the direct transmission matrix. The identification step is quite similar, 
if not identical, to the identification of a linear system [1,2]. This step establishes a 
specific set of coordinates for the whole identification process. This set of coordinates is 
not unique, depending mainly on the size of Hankel matrix and the resulting choice of 
matrix that represents the observability matrix. The second step starts by generating 
another set of pulse responses with the same constant input as the frst step but for 
multiple sample periods. These multiple-pulse responses are used to define another set of 
Markov parameters to form a Hankel- like matrix for each input. The observability matrix 
obtained in the first step is then applied to the Hankel- like matrix to compute the 
corresponding controllability matrix of the input to identify the input vector and the 
coefficient matrix associated with the coupling terms between the state and the input. 

Simple examples are given to demonstrate how to apply the method to identify a 
continuous-time bilinear system and how to transfer the identified model from one set of 
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coordinates to the other set of coordinates. The coordinate transformation also serves as 
a way of verifying the identified system. 

Basic Formulation 

Let x and A c be the state vector of dimension n x 1 and its corresponding state 
matrix of n x n , and u and B c be the input vector of r x 1 and its corresponding input 
matrix of nxr . Subscript c signifies the associated quantity in the continuous-time 
domain. The bilinear state equation in the continuous-time domain is commonly 
expressed by 

x = A c x+ B c u+ Y N d xu t ( 1 ) 

<=i 

where the coupling term xu, between the state vector x and each individual w, (/= 1,..., r) 
in the input vector u is weighted by the matrix N a of n x n . The measurement equation is 
identical to the one for a linear system that is commonly described by 

y=Cx+Du (2) 


where y is the output measurement vector of m x 1 , C is the output matrix of mxn and D 
is the direct transmission matrix of mxr . 

For simplicity, consider only one input at a time. Equation (1 ) reduces to 


x = A c x +b ci u j + N d xu i 


(3) 


where b c \ is the /th column of B c associated with the input u^. Assuming u t =v> i where n, 
is a pre-specified constant, the continuous-time state equation (3) further reduces to 


x = (A c +v i N.,)x + b a v i 
The discrete-time model of this system is 

x(k + 1) = A l x(k) + b i : i = l,2,...,r 
with the measurement equation 


y l (k) = Cx(k) + d i 

where d x is the /th column of D associated with the input tik and 

__ N n ) A/ 





dx b c p, 


(5) 

( 6 ) 

(7) 

( 8 ) 
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d i = d^) i (9) 

The quantity A t is the time interval for data sampling. Assuming that the initial state 
jc( 0) is a zero vector, i.e., x(0) = 0„ xl , the state response for the discrete-time model, Eq. 


(5), can be computed by: 

*( 0 ) = 0 HX1 

x( \) = A i x(0) + b i =b i 

x(2) = A 2 x(\)+b i = A i b i + b i 

x(3) = A?x(2) + E = A;b. + Ah + b k 

x(N) = A j N x(l)+b i =A i N '% + A 1 N ~% J f-+Ap i +b f 


After / > NAt , let »,(/) = 0. The state equation (3) reduces to the simple form 


II 

(11) 

Its discrete-time model is 


x(k + 1) = Ax(k) 

(12) 

where 


<5 

II 

(13) 

The free decay response after t > NAt becomes 

x(N + 1) = Ax(k) = A (A^'% + A, N ~ 2 b t + • • • + Aft +Z>.) 
x(N + 2) = A 2 x(k) = A 2 ( 4 N ~'b, + 4"~ 2 E +--- + A,b i +b i ) 

(14) 

x(N + {) = A e x(k) = A c ffl- , li.+ 4"- 2 bi+'-- + 4b I +bj) 



where C is an integer indicating the data length of the free-decay response. 

From Eqs. (2), (10), and (14), the measurement quantities y,(k) for 

k = 0,1, • • •, N + C due to the force excitation of u k = o A can thus be computed as 
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v,( 0) = Cx(0) + d i =d l 

y,(l) = Cx(l) + 4 = Cb, +d f 

v,(2) = Cx(2) + d. = C(Afy + b i J+d i 

y i (N) = Cx(N)=C(A i N -%+A, N - 2 b i +---+Afr+b i ) ° 5) 

y , ,(N + 1 ) = Cx( N + 1) = CA + A^'b, +--- + Afr+b l ) 

y, ( N + C) = Cx(N + C) = CA'( f~ l b, + A^b, +--- + Afy+b i ) 

The upper portion, v(0), **, y(N) , of Eq. (15), corresponds to the multiple-pulse 
response resulting from a constant force over multiple sample periods. But the lower 
portion, y(N + l),y(N + 2),---,y(N + (), corresponds to the free-decay response. 

It is clear that the system matrices/vectors, A c , b d , N ci ,C,d j , are embedded in the 
output quantities shown in Eq. (15). First of all, we need to use the output measurements 

to extract the discrete-time matrices/vectors, A, C, 4 , b n d r It is worth to stress that the 

multiple -pulse response and the tree-decay response result from two different discrete 
models. 

The free-decay response, y(N + 1 ), v( N + 2),- • -,y(N + C) , after t > NAt is quite 
similar, if not identical, to the pulse response for a linear system. Any linear system 
identification technique may be applied to compute the state matrix A and the output 
matrix C. The key idea is to make judicious use of this linear portion of the bilinear 
system. The identification problems for linear systems have been extensively studied and 
many good techniques have been developed and implemented. 

System Identification Method 

The identification method requires two steps. The first step is to identify the state 
matrix A c , the output matrix C, and the data transmission matrix D. The second step is to 
determine the input matrices B c , and N\ for the coupling term between the state vector x 
and the ith input u\. 
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Identification ofA c , C, and D: 


First, let us apply a pulse of magnitude i), to the system for one time step At to 

generate the pulse response for the first input u\. From Eq. (15) for N = 1, the pulse 
response has the following expression. 

v,(0) = d l 

y l (\) = Cb ] 

y,(2) = CAb l (16) 


y,(f + 1) = CAfy 

All other input pulse responses can be similarly generated to yield 

v 2 (0) = <7 2 ••• v,.(0) -d r 

v 2 (l) = Cb 2 ••• v,.(l) = Cb r 

v 2 (2) = CAb 2 ••• y,.(2) = CAb r (17) 

y 2 (C + 1) = CA r b 2 ••• y r (C + \) = CA ( b r 

Equation (17) is obtained by replacing the subscript 1 representing the first input with the 
other input integers 2 through r. Let us define the system Markov parameters to be 

^I(O) = [ v,(0) y 2 ( 0) - j r (0)] = [j, d 2 - d r ] = D 

W = b,(D y 2 ( i) - yM=[cb, cb 2 - cF] 

^(2) = [y,(2) y 2 ( 2) - y,.(2)] = [CAb i CAb 2 - CAb r ] (18) 


})(f+l) = [y,(f + 1) y 2 (f+ 1) ••• v,(( + l)] = [c4 CA r b 2 ••• C4b r ] 

The use of the subscript 1 for Y\(k) ( k = 1,2,..., C + 1 ) is intended to signify one-time-step 
pulse response. Equation (18) provides the basic parameters for system identification. 
Indeed, let us form a Hankel matrix as follows. 


ff,= 

'^(1) 

Y](2) 

YX 2) - 

^(3) 

W) 

- YX\ 3+1) 

= 

C 

CA 

[«, AB, ■ 

•• A^B,]{\9) 


JXCL) 

>^(a+l) •• 

•• Wa + P-1)_ 


CA a ~\ 




where 
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£,=[>, b 2 ••• b] (20) 

The matrix product on the right-hand side of Eq. (19) shows the relationship between the 
system Markov parameters and the discrete-time system matrices. Obviously the Hankel 
matrix H\ has the rank n that is the order of the state matrix A if we choose a and P 
such that a wand Pr are larger than or equal to n where m is the number of outputs and r 
is the number of inputs. Using the singular value decomposition (SVD) to decompose 
the Hankel matrix H\ yields 

( 21 ) 


where Z,of nxn is a square matrix containing n non-zero singular values. The matrix 
(7, is of dimension amxn and the matrix V x is of dimension P rxn. 

From Eq. (19), one may choose 


and 


C 

CA 


U t = 


CA a ~ l 


( 22 ) 


AB x ••• A*-'B X ] (23) 

This choice is not unique. Many other choices are also valid. The other common choice 
is 


uX !2 


C 

CA 


CA a ~' 


(24) 


and 


' 2 V X T = [% AB X ••• (25) 

Note that the choice of Eq. (22) has the advantage that 

u!u t =l,„^u!=u; (26) 

because (/, is a unitary matrix resulting from the property of singular value 
decomposition. Nevertheless, the choice of Eqs. (24) and (25) has a nice property of 


7 



balanced coordinates. Equation (22) or (24) is commonly called observability matrix 
whereas Eq. (23) or (25) is referred to as the controllability matrix. 

Equations (22) and (23) produce the following solutions 

C = the first m rows of t/, (27) 

B t = the first r columns of S,F, r (28) . 

Since the choices of controllability and observability matrices are not unique, the 
identified matrices C and B ] are not unique. To determine the state matrix A, let us first 
define and observe the following matrices. 

C 

CA 

= . (29) 

CA"~ 2 _ 

and 


= U^A (30) 

Deleting the last m rows of (/, forms the matrix t/, T whereas deleting the first m rows of 
(/, yields the matrix U a . It is then clear that the state matrix A can be determined by 

A = U' it U a (31) 

For the identified state matrix to have the rank n , both (a -1 )m x n matrices t/ |T and U ii 
must also have the rank n. This implies that a must be chosen such that (a -1 )m > n, 
i.e., am > n . Of course, we have assumed that the pulse force u, for i = 1 , 2 ,..., r are 
chosen so that all system modes are excitable and observable. 

With the aid of Eq. (13), Eq. (31) produces the continuous -time state matrix as 

^ilogt/O^loglt/X) (32) 




CA 

CA 2 

ca r 1 


Note that the conversion from a discrete-time state matrix to a continuous-time state 
matrix may not unique. 


i 
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To this end, we have determined A c from Eq. (32), C from Eq. (27), 5, from Eq. 

(28), and D from Eq. (18). The original transmission matrix D can be recovered using 
Eq. (9) to have 

D=D diag[l/i), 0> 2 - 1 /d,.] (33) 

Let us stress that the identified matrices A c , B i and C are not uniquely determined but D 
is coordinate invariant and so is uniquely determined. 

Identification of B c , and N ci : 

The second step begins with generating the two-sample-period pulse response for 
all inputs with one input at a time, i.e., a force is applied with the same magnitude as 
above to the system for two time steps 2 At . From Eq. ( 1 5) for TV = 2, we obtain 

v,(0) = d x ••• y r (0) = c/ r 

y^^Cb.+d, - v,.(l) = Cb r +d r 

Ti( 2) = C[^ + ^] ••• y,.(2) = C[4^. +b r ~\ (34) 

y l {C + 2) = CA ( [A l b i +b i '] ••• y r (C + 2) = CA! \lj> r + b r ] 

Now define the system Markov parameters from the two-sample-period pulse response as 
T 2 (2) = [v,(2) *(2) - v r (2)] 

= [c(^+ij) C(l 2 b 2 +b 2 ) - c(M+F r )] 

n(3) = b,(3) y 2 ( 3) y,.(3)] 

= [CA(Ah+b) CA(A 2 b 2 +b 2 ) ••• CA{A,b r + b,)\ (35) 
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where 


B 2 =[{A,b,+b x ) (A 2 b 2 + b 2 ) ••• (4A +/>,.)] (37) 

With the help of Eq. (22), the nxr matrix B 2 in Eq. (36) can be solved by 


b 2 = 


c 

CA 


CA 


H-l 


h 2 =u;h 2 


(38) 


Similarly, we may continue the process to generate three-sample-period pulse 
response, four- sample -period pulse response, etc. up to the p-sample-period pulse 
response for all inputs with one input at a time using a force of the same magnitude as 
earlier applied to the system for p time periods p At . From Eq. (15) for N=p, we have 
y l (p)=C(lr%+- + b ] ) ••• y r (p) = C(A, p ~ l b r +--- + b r ) 

y i (p + \)=CA(Ar%+- + b l ) ••• y r (p + l) = CA(A 1 p -'b r + - + b r ) 
y x {p + 2)=CA 2 {A p ~%+-+b x ) ••• y r (p + 2) = CA 2 (A! , -%+"' + b r ) (39) 

y ] (p-K) = CA , (AT%+- + b i ) - y r {p + C) = CA ( (AT'b r +-+b r ) 


Now define the system Markov parameters for the p- sample- period pulse response as 

Y P (p)=[y\(p) y 2 (p) — y r {p)\ 

— x b x +---+b l 'j C^Af'b-,-1 h b 2 j C(A pl b r - 1 

y p (p+i) = ly t (p+i) y 2 (p+ 1) - y r (p+ 1)] 

= \ L CA(A p -% + ~- + b x ) CA[A p -%+--- + b 2 ) ••• CA(A, pA b r + ••• + £)] ( 40 ) 

y p (p+0 = [y(p+0 y 2 (P+ 0 - y r (P+0] 

= [CA'(A x p -%+-+b l ) CA' (A p ~'b 2 +--- + b 2 ) ••• CA f + --- + b r )] 


Let us form a amxr matrix as follows. 



y P (p) 


C 

H »= 

Y p (P + () 

= 

CA 


Y p (P + a-1) 


CA n ~'_ 


( 41 ) 
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where 


B ,= 


_(Ar'b l +-+b l ) (Af'b 2 +---+b 2 ) ••• [A? % +•••+/>,.) 
With the help of Eq. (22), the n x r matrix B p in Eq. (41 ) can be solved by 


(42) 


B P =\ 


C 

CA 


H,=U\H p 


(43) 


CA 


(X — 1 


To determine B c , let us first observe the matrices B { , B 2 ,---,B p defined in Eqs. 
(20), (37), and (42), and determined by Eqs. (28), (38), and (43), i.e.. 


5=[ 

/?, b 2 ••• b r ] 


B 2 = 

(Afy+b,) (l 2 b 2 +b 2 ) ••• (A r b r +/T)] 

(44) 


\Ar%+-+b, ) (A?-%+- + h 2 ) ■■■ (Ar%+-+b,)\ 


Applying the recursive formula 




(45) 

yields the controllability- like matrices, 



C,=[6, Ah = 

(46) 

To determine the state matrix A j , let us first define the two matrices 



C^=[», 4A - -r 2 *,] 

(47) 

and 

C„-[^ ?A 

(48) 


Deleting the last r columns of C, forms the matrix C (< _ whereas deleting the first r 
columns of C, yields the matrix . Equations (47) and (48) produce the solutions: 


b t = the first r columns of C, (49) 

A, = C,_C^ (50) 
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for i = 1, 2,..., r. For the identified matrix A f to have the rank «, both nx(p- 1) 
matrices C (< _ and C must also have the rank n. It implies that p must be chosen such 
that p — 1 > n . This indicates that the system identification method requires a total of at 
least (n +1 ) sets of responses generated by (« +1 ) different time periods of pulse input. 

Based on Eqs. (7) and (8) for the definitions of A { and b t , taking the conversion 
from discrete- time to continuous -time produces 

4 +»,v, ( =-l|og(4) = i|og(C,jCL) (51) 

and 


h a = — U +tt(4 +t ) A,)(A/) 2 +— (A c +v> i N a ) 2 (Atf +■■■ h 


for i = 1 , 2 ,. . r where I nX(l is a nxn identity matrix, that, in turns, yields 

4 =(4i ^<2 •" K ] (53) 

Again, one should be cautious to take the conversion because of its non- uniqueness 
problem [1], From Eqs. (32) and (51), the matrices N d ( i= 1,2,..., r ) are determined by 




To this end, we have identified all continuous-time system matrices A c , B c , N c - h C, 
and D for the bilinear system described by Eqs. (1) and (2) from pulse responses 
generated by pre-specified pulse inputs. 


Coordinate Transformation 

Let the state vector x of nx 1 in Eqs. ( 1 ) and (2 ) be transformed to the new state 
vector x of «xl by the nonsingular transformation matrix <t> of nxn. Equations (1) 
and (2) become 

x = A c x+ B c u + X N ci xu j (55) 

/=i 

and 

y= Cx+ Du (56) 

where 
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( 57 ) 


x = <t>-'x 
A. =<D"'4.<D 
N d =<P'N a <\> 

B = 0 > 'B 

c c 

c = c<t> 

It is clear that the transformed matrix A c is similar to the original matrix A c in the sense 
that their eigenvalues are identical. The same statement is also true for the matrices N d 


and N d . One question that may arise is, for given two sets of matrices representing the 

same bilinear system, what the transformation matrix should be to convert from one 
coordinate to the other. 

First, form the observability matrices for both sets of system matrices 


Q = 


C 

CA. 


CA 


«-i 


(58) 


and 


Q = 


c 

CA„ 


«-! 


CA 


Substituting the relationship from Eq. (57) yields 


(59) 



c 


c 

Q= 

CA c 

= 

CA 

C 


ca 


% 

1 


<J> = 0O 


(60) 


that in turn provide the following solution for computing the transformation matrix <t> 


® = Q f Q 


(61) 



This transformation matrix <J> will transform the original system coordinate to another 
system coordinate without changing the input-to-output map. 


Numerical Example 


Consider the following example presented in Ref. [ 4 ] 


where 


x = A c x+ B c u + N cl xu l + N c2 xu 2 
y = Cx 



"-1 0 " 


O 

o 


'i r 

A = 


; N, = 


; AC, = 


c 

1 -2 

’ ti 

1 1 

’ c 2 

0 0 


B = 


1 0 
0 1 


; C = [0 1] 


(62) 


(63) 


Assume that we don’t know the order of the system. Let us generate five sets of data 
with the time period A t = 1 second. The first set of data with one input at a time is 
generated by a unit force of period 1 second. The second set of data with one input at a 
time is obtained by applying a unit force of period 2 seconds and the fifth set of data is 
computed with a unit force of period 5 seconds. The total of data points is set to be 20 
for each data record. 

Figure 1 shows a total of 10 ( p xm xr ) responses from 2 inputs (r = 2), 1 output 


(m =1), and 5 different multiple-pulse inputs (p = 5). Each response sampled at 1 Hz has 
20 data points. These ten responses are obtained by numerically integrating the bilinear 
system shown in Eq. (62). With a = 5 and 0 = 6 , the Hankel matrix H\ shown in Eq. 
(19) should have the size of 5x12 (amx(3r). The state matrix and the output matrix 


identified from this Hankel matrix is 


A.. = 


-1.0629 

0.0148 


3.9782 
-1.9371 ’ 


C = [-0.9355 0.3497] 


(64) 


The singular values of this Hankel matrix are 

=diag [0.8347 0.0543 0 0 0] 

implying that the order of the system is n = 2. The other matrices Hy for k = 2, 3, ...,p 
shown in Eqs. (36) and (41 ) are of the size of 5x2, that produce the matrices 
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5, , B 2 ,•••, B 5 of 2 x 2 each shown in Eq. (44), and in turn yield C, and C 2 of 2 x5 
each shown in Eq. (47). The B c , N c \, and N C 2 can thus be identified as 
- f 1.7752 3.29111 - f 0.1678 0.31111 


1.7752 

3.2911" 

N = 


l ' [-0.4182 

-0.7752_ 

- [ -0.0929 

-0.9824" 

B c = 


1 -0.2484 

0.2314 


; K-2 = 


0.4489 0.8322 


The quality of the identification of the system is evaluated next. The following 
transformation matrix is computed from Eq. (6 1 ) 

-0.8715 -3.6998" 

0 = ( 6 < 
-0.9355 0.3497 

the matrices A., B c , C , N c] and N c2 would be transformed by using Eq. (57) to become 

A c , B c , C, N c \, and N C 2 shown in Eq. (63). 

Let us examine another case where we keep the unit force for the first input but 
change to 0.5 for the second input. Figure 2 shows the 10 multiple-pulse responses. 
Applying the same identification procedure as above, we obtain the following system 
matrices 


-0.9509 3.9503 

0.0131 -2.094 


C= [-0.9253 0.3760] 


j. _ 1.8543 3.21681 ~ [0.1898 0.3292" 

~~ |_ -0.4925 -0.8543J’ ‘-"[0.4670 0.8102 

- [-0.0998 -0.97551 

B ;.= 

‘ [-0.2457 0.2591 

with the transformation matrix 

[-0.9757 -3.67371 

0 = 

-0.9253 0.3760 

Note that the Hankel singular values are 

I, =diag[0.6753 0.0489 0 0 0] 

This set of system matrices also represents the bilinear system, because it can be 
transformed back to the original system exactly. 
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Concluding Remarks 


A new method is introduced for identification of a continuous-time rrulti- input 
and multi-output bilinear system. The approach is to make judicious use of the linear- 
model properties of the bilinear system when subjected to a constant input. It has been 
shown in this paper that a bilinear system can be treated as a combination of two linear 
systems in the identification process. The first linear system is the one obtained by 
deleting the nonlinear terms of the bilinear system. The second linear system is given by 
assuming a constant input. Due to this latter property, the identification process for the 
bilinear system becomes a combination of two linear-system identification processes. 
The key is to combine these two linear-system identification processes in the same 
coordinate system. The resulting identified system matrices would be similar to the 
original ones in the sense that they represent the same bilinear system but in different 
coordinates. With a proper coordinate transformation, both the original model and the 
identified model are identical. 
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